#include "auxiliar.h"

void printOnFile_Cij(ofstream &fid, pMesh theMesh, const int &meshDim, set<int> &setOfDomains);
void printOnFile_NodalVolume(ofstream &fid, pMesh theMesh, const int &meshDim, set<int> &setOfDomains);
void printOnFile_Dij(ofstream &fid, pMesh theMesh, const int &meshDim, set<int> &setOfDomains);
bool getDataAssociatedToEntity(pEntity ent, int dim, int dom, char** str, double *buffer);
void printfData(ofstream &fid, pEntity ent, int dim, const double *buffer, char** str); // char** str -> BACALHO!! adicionado em todas
                                                                                        // as funções de impressão para poder saber
                                                                                        // se é coeff. C ou D!! MUDAR ISTO DEPOIS!
                                                                                        // BRUNO LUNA (29/05/10)

// the following function exports geometric coefficients based on finite volume
// method described on Darlan thesis only. Coefficients are: Cij , Dij and nodal
// volume
// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
// IN
// theMesh - pointer to mesh object
// meshDim - mesh dimension 2D/3D
// setOfDomains - set container containing flags associaated to all sub-domains
// OUT
void exportCoefficients(pMesh theMesh, const int &meshDim, set<int> &setOfDomains, const char *meshFilename)
{
	cout << "export coefficients\n";

	ofstream fid;
	fid.open(meshFilename);
	fid << setprecision(18) << fixed;

	// print in the begining of file all flags associated to all domains
	// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
	set<int>::iterator iter;
	fid << setOfDomains.size() << endl;
	for (iter=setOfDomains.begin(); iter!=setOfDomains.end(); iter++) fid <<*iter<<" ";
	fid << endl;

	printOnFile_Cij(fid,theMesh,meshDim,setOfDomains);
	printOnFile_NodalVolume(fid,theMesh,meshDim,setOfDomains);
	printOnFile_Dij(fid,theMesh,meshDim,setOfDomains);

	fid.close();
}

// printf on file Cij coefficients
// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
void printOnFile_Cij(ofstream &fid, pMesh theMesh, const int &meshDim, set<int> &setOfDomains)
{
	double buffer[3];
	char  *cdstr[3] = {"Cx-","Cy-","Cz-"};

	pEntity edge;
	set<int>::iterator iter;
	// list Cij for domains
	for (iter = setOfDomains.begin(); iter != setOfDomains.end(); iter++)
	{
		fid << "domain:"<<*iter<<endl;
		int edgeID = 0;
		EIter eit = M_edgeIter(theMesh);
		while (edge = EIter_next(eit))
		{
			edgeID++;
			// check if Cij exist for domain=(*iter). if true, print it on file
			if (getDataAssociatedToEntity(edge,meshDim,*iter,cdstr,buffer))
			{
				fid << edgeID <<" ";
				printfData(fid,edge,meshDim,buffer,cdstr);       // cdstr -> BACALHO!! adicionado em todas
                                                                 // as funções de impressão para poder saber
                                                                 // se é coeff. C ou D!! MUDAR ISTO DEPOIS!
                                                                 // BRUNO LUNA (29/05/10)
			}
		}
		EIter_delete(eit);
		fid << "end\n";
	}
}

// printf on file volume of control volume (nodal)
// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
void printOnFile_NodalVolume(ofstream &fid, pMesh theMesh, const int &meshDim, set<int> &setOfDomains)
{
	set<int>::iterator iter = setOfDomains.begin();
	for (; iter != setOfDomains.end(); iter++)
	{
		fid << "domain:"<<*iter<<endl;
		VIter vit = M_vertexIter(theMesh);
		while (pEntity vertex = VIter_next(vit))
		{
			// total nodal volume
			double val = .0;
			char volumeString[64];
			sprintf(volumeString,"volume-%d",*iter);
			EN_getDataDbl(vertex,MD_lookupMeshDataId(volumeString),&val);

			if (val > .0)
			{
				fid  << EN_id(vertex) << " "<< val << " ";

				// weighted porosity
				val = .0;
				EN_getDataDbl(vertex,MD_lookupMeshDataId("wght-porosity"),&val);

				fid  << val << endl;
			}
//			if (getDataAssociatedToEntity(vertex,meshDim,*iter,vertexStr,buffer))
//				printfData(fid,vertex,meshDim,buffer);
		}
		VIter_delete(vit);
		fid << "end\n";
	}
}

// printf on file Dij coefficients: associated to boundary edges
// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
void printOnFile_Dij(ofstream &fid, pMesh theMesh, const int &meshDim, set<int> &setOfDomains)
{
    // double buffer[3];
    //char  *cdstr[3] = {"Dx-", "Dy-", "Dz-"};
    // #ifdef MULTIGRID
    double buffer[6];
    char  *cdstr[6] = {"Dx-", "Dy-", "Dz-","deltaDx-","deltaDy-","deltaDz-"};
    // #endif
	set<int>::iterator iter;

	if (meshDim == 2)
	{
		pEntity edge;
		for (iter = setOfDomains.begin(); iter != setOfDomains.end(); iter++)
		{
			fid << "domain:"<<*iter<<endl;
			EIter eit = M_edgeIter(theMesh);
			while (edge = EIter_next(eit))
			{
				if (getDataAssociatedToEntity(edge,meshDim,*iter,cdstr,buffer))
					printfData(fid,edge,meshDim,buffer,cdstr);                          // cdstr -> BACALHO!! adicionado em todas
                                                                                        // as funções de impressão para poder saber
                                                                                        // se é coeff. C ou D!! MUDAR ISTO DEPOIS!
                                                                                        // BRUNO LUNA (29/05/10)
			}
			EIter_delete(eit);
			fid << "end\n";
		}
	}
	else
	{
		// printf on file Dij coefficients: associated to boundary faces
		// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
		pEntity face;
		pGEntity gEntity;
		for (iter = setOfDomains.begin(); iter != setOfDomains.end(); iter++)
		{
			//fid << "domain:"<<*iter<<endl;
			FIter fit = M_faceIter(theMesh);
			while (face = FIter_next(fit))
			{
				gEntity = (pGEntity)F_whatIn(face);
				int faceFlag = GEN_tag(gEntity);
				if (getDataAssociatedToEntity(face,meshDim,faceFlag,cdstr,buffer))
				{
					printfData(fid,face,meshDim,buffer,cdstr);                          // cdstr -> BACALHO!! adicionado em todas
                                                                                        // as funções de impressão para poder saber
                                                                                        // se é coeff. C ou D!! MUDAR ISTO DEPOIS!
                                                                                        // BRUNO LUNA (29/05/10)

				}
			}
			FIter_delete(fit);
			fid << "end\n";
		}
	}
}

bool getDataAssociatedToEntity(pEntity ent, int dim, int dom, char** cdstr, double *buffer)
{
	int i;
	double sum = .0;

	dim = ( ent->getType()==mEntity::VERTEX )?1:dim;

	//char *strbuffer[dim];
	//for (i=0; i<dim; i++) strbuffer[i] = new char[64];

    //#ifdef MULTIGRID
    char *strbuffer[6];
	for (i=0; i<6; i++) strbuffer[i] = new char[64];
    //endif MULTIGRID

	for (i=0; i<dim; i++)
	{
		buffer[i] = .0;
		sprintf(strbuffer[i],"%s%d",cdstr[i],dom);
		EN_getDataDbl(ent,MD_lookupMeshDataId(strbuffer[i]),&buffer[i]);
		sum += fabs(buffer[i]);
	}

	if(cdstr[0]=="Dx-")
    {
        //#ifdef MULTIGRID
        for (i=3; i<(dim+3); i++)
        {
            buffer[i] = .0;
            sprintf(strbuffer[i],"%s%d",cdstr[i],dom);
            EN_getDataDbl(ent,MD_lookupMeshDataId(strbuffer[i]),&buffer[i]);
        }
        //endif
    }
	return (fabs(sum) > 1e-15);
}

void printfData(ofstream &fid, pEntity ent, int dim, const double *buffer, char** cdstr)
{
	if ( ent->getType()==mEntity::VERTEX )
		fid <<EN_id(ent) << " " << buffer[0] << endl;
	else
	{
		int i;
		vector<pVertex> node;
		M_GetVertices(ent,node);
		for (i=0; i<node.size(); i++) fid << EN_id(node[i]) << " ";
		node.clear();
		for (i=0; i<dim; i++) fid << buffer[i] << " ";

		if(cdstr[0]=="Dx-")
        {
		// #ifdef MULTIGRID
		for (i=3; i<(dim+3); i++) fid << buffer[i] << " ";
		// endif
        }
		fid << endl;
	}
}
